Theoretical Population Biology

Price equation and genotype dynamics are two methods for studying the fixation of one allele by natural selection in a diploid population. There are two strict monotonicity conditions that imply the fixation of one allele. The genotype dynamics is called Haldane monotone if the relative frequency of one allele strictly increases along all solutions of the genotype dynamics, so this allele is fixed. In this paper, we show that the genotype dynamics is Haldane monotone if and only if the right-hand side of the Price equation is always strictly positive. The other strict monotonicity condition requires that the relative frequency of a homozygote strictly increase according to the genotype dynamics. For example, in a model where the genotype dynamics is governed by interactions between individuals, the cost-accepting homozygote is fixed by natural selection if the other genotypes always receive a smaller average gain from all interactions than the cost-accepting homozygote. Both monotonicity conditions require that the interaction is not well-mixed in the population. These two conditions are not equivalent. In addition, we give a non-monotonicity condition, which also implies the fixation of a homozygote. The fixation of a homozygote depends on the phenotypic payoff of the interaction, the genotype-phenotype mapping, and the interaction scheme. In a sexual population, the interaction scheme of siblings depends on the mating system, and so do the conditions of fixation of the cost-accepting homozygote. We present examples showing that if we only change the monogamous mating system, assuming panmixing or mating assortativity, then the condition for the fixation of the cooperator homozygote is b > 2c and b > c, respectively.

Interest in quantifying linkage disequilibrium (LD, non-random associations of alleles at different loci) has skyrocketed in recent years as researchers have focused on use of LD in genome-wide association studies (GWAS), for studying historical demography, and for estimating effective population size (Ne). The most widely used LD metric is r2 = the squared correlation of alleles at a pair of loci. Despite a half century of efforts, developing an unbiased expectation of r2 as a function of the many factors that can affect it (physical linkage, genetic drift, selection, migration, mutation, mating systems) remains elusive. Furthermore, even when all of these other factors are absent, empirical estimates of r2 are upwardly biased by sampling a finite number (S) of individuals, and that must be accounted for if one wants to focus on the desired signal of LD. Previous approaches to estimate [Formula] have been shown to be biased to greater or lesser degrees. The purpose of this short paper is to demonstrate that a simple and apparently exact expression for [Formula] does exist for the special case where sampling error is the only factor contributing to r2, in which case [Formula] = 1/(S - 1). When other factors contribute heavily to LD, [Formula] shrinks toward 0 as empirical r2 [->] 1. However, for estimating contemporary Ne with unlinked markers, empirical r2 will generally be small and 1/(S - 1) will provide a robust estimate of [Formula].

A sudden reduction in population size increases the rate of genetic drift, reducing variability and increasing the mean level of homozygosity. The resulting increased exposure of recessive or partially recessive, strongly deleterious alleles to selection against homozygotes may lead to their being purged from the population, potentially allowing mean fitness to increase after an initial decline, and accelerating the decline in inbreeding depression associated with reduced variability. However, detailed population genetic theory on the effects of population bottlenecks on mean fitness and inbreeding depression remains limited. We develop a theoretical framework for small, randomly mating populations founded from a large population near mutation-selection-drift equilibrium, using both simulations and approximate analytical predictions. These provide quantitative predictions for the dynamics of the populations mean fitness and level of inbreeding depression following a bottleneck. In particular, we derive an approximate expression for the time needed for mean fitness to recover after an initial decline; such a recovery requires selection to be sufficiently strong relative to drift and mutations to be sufficiently recessive. In contrast, weakly deleterious mutations cause reductions in mean fitness and inbreeding depression that are similar in size to those predicted from increases in neutral homozygosity.

Darwinian fitness is equated here with invasion fitness and defined as the quantity determining the fate--certain extinction or possible spread--of a single mutant type. We derive it, together with its phenotypic derivative, for evolution in group-structured populations under limited genetic mixing, where the demography of the focal species and its environment is modeled as a discrete-time stochastic process. Reproduction, physiological development, dispersal, and survival are influenced by interactions within and between groups and by environmental fluctuations within and across generations. Using multitype branching processes in random environments, we show that invasion fitness is predicted by a stochastic growth rate that can be represented biologically in two meaningful genealogical ways. First, as the long-term geometric mean of the expected per-capita number of mutant copies produced per time step by a representative member of the mutant lineage. Second, as the the long-term geometric mean of the expected reproductive-value-weighted per-capita number of mutant copies produced by such an individual. This latter representation is useful for computing the phenotypic directional derivative of invasion fitness. Moreover, this derivative can be written as an actor-centered inclusive-fitness effect derived from properties of the resident population process. This effect depends on class-specific fitness differentials, relatedness, reproductive values, and class frequencies. However, unless generation- and class-specific fitness defines a stochastic matrix, the derivative does not separate stochastic reproductive values from relatedness and class frequencies, and must be evaluated by simulations. In summary, we formalize invasion fitness biologically quite generally and show how Hamiltons marginal rule is deduced from it.

Conventional LD measures such as r2 perform poorly in the rare common regime, particularly in asymmetric configurations such as nested haplotype structure. Because r2 is symmetric and quadratic, it removes directional structure in two ways: squaring discards the sign, or phase, retained by the signed LD coefficient D, while symmetric normalization hides the asymmetry between the conditional probabilities P(A|B) and P(B|A). Although D recovers the phase, it is locus symmetric and unnormalized; its magnitude is hard to compare across frequency regimes and it does not by itself express which way the asymmetry runs. We therefore analyze the conditional-probability asymmetry {Delta} = P(A|B) - P(B|A), together with r2 and D, as distinct scalar functions on the haplotype simplex under the Fisher information metric. The conditional probabilities P(A|B) and P(B|A) are bounded in [0, 1], directly express carrier-set inclusion, and are more readily visualized than D. Moreover, their difference admits the exact decomposition {Delta} = M + C into a marginal frequency term M and an LD-coupled term C. Prior work has characterized either the mathematical behavior of LD normalizations across allele-frequency space or the Fisher geometry of the haplotype simplex, but not their connection. We bridge this gap by showing that the geometric structure of the simplex explains why LD measures disagree in the rare common regime and why symmetric normalizations such as r2 lose directional information. We show that the fixed-frequency leaf is intrinsically anisotropic, positively curved, and frequency-dependent under the Fisher metric. These geometric predictions are tested empirically , in phased 1000 Genomes data1 and a two locus Wright Fisher model, in a companion paper (Ichikawa, preprint); the present note develops the geometry itself. Keywords: linkage disequilibrium; Fisher information metric; haplotype simplex; rare variant; conditional-probability asymmetry; nested haplotype structure

Despite sharing the same genes and the same environment, individuals often develop substantial phenotypic differences. While this pattern has been documented across diverse species and traits, the processes giving rise to this "stochastic" or non-shared environmental variation remain unclear. Recent mathematical models of development in which phenotypes are gradually constructed may offer some clues. These models show that imperfect environmental cues can generate striking variation in developmental trajectories and adult phenotypes. At the population level, such imperfect cues produce increasing stability of individual differences across ontogeny (e.g. animal personality) and patterned distributions of mature phenotypes (e.g. normal or skewed) that resemble those observed in real organisms. Our paper synthesizes existing models in which stochastic phenotypic variation arises solely as a by-product of mechanisms missing their phenotypic targets because of imperfect cues. We then link these models to related, but independent, mathematical theory exploring the environmental conditions under which stochastic phenotypic variation is favoured by natural selection. Our integration shows that stochastic sampling is often favoured over classic bet-hedging strategies involving non-plastic generalist or specialist strategies. Our findings provide new directions of research on stochastic sampling as a mechanism for adaptive stochastic variation within and across generations.

1Although resources are typically distributed continuously in space, species distributions often organize into discrete clusters. In his seminal paper [36], Turing demonstrated that such clusters can spontaneously arise in population densities, even when populations evolve in environments with continuously varying conditions. This phenomenon is known as Turing instability. In this work, we focus on two models grounded in population dynamics: a one-dimensional model based on the nonlocal Fisher-KPP equation, and a two-dimensional model involving an environmental gradient. We show that phenotypic clusters (sometimes referred to as "species") emerge in these models. We prove that they do not emerge because of Turing instability, but because of stochasticity, and that they disappear when stochasticity is reduced. First, for both models, we start our simulations with initial populations uniformly distributed in the state space. We show that phenotypic clusters quickly emerge and that the distances between them depend on the population size, that is, on the degree of stochasticity. Next, we start from already clearly defined phenotypic clusters. We identify three regimes in the connection between population size, the initial distances between clusters, and the distances between clusters at equilibrium. Last, on the two-dimensional model, we relax the hypothesis of complete clonality by varying the effective recombination rate, explore its effect on phenotypic clustering, and show that phenotypic clustering decays drastically with slight recombination.

A phylogenetic network represents evolutionary relationships involving hybridization, gene flow, or admixture. While the full network may not be identifiable from genomic data under common coalescent models, its tree of blobs, depicting only the tree-like portions of the network structure, is. We introduce ECToBlob (Edge Contraction for Tree of Blobs), a new statistically-consistent algorithm to estimate the tree of blobs from quartet concordance factors. Starting from a resolved tree, ECToBlob successively contracts edges which statistical tests indicate do not belong in the tree of blobs, due to reticulate or polytomous signal. We show that ASTRAL provides a valid starting tree under common assumptions, in that, asymptotically in the number of loci, trees optimizing ASTRALs criterion refine the tree of blobs. We describe several algorithm variants, differing in how evidence from multiple tests are combined to determine if the edge should be contracted, and provide software implementations. Relevance to Life SciencesHybridization, gene flow, or admixture are now recognized as important aspects of evolutionary history, but their genomic signal is confounded with that from a coalescent process, creating substantial challenges for inferring phylogenetic networks. The networks tree of blobs identifies areas where reticulation occurred, separated by tree-like branching. ECToBlob quickly estimates the tree of blobs using quartet concordance factors from gene trees, and provides a measure of statistical support for its result. Performance is illustrated through simulation and on empirical data, using an implementation in the R package MSCquartets. While the presence of a blob may be all that can be inferred in some cases, in others ECToBlob offers a robust and principled way to focus further analyses on more local reticulate structure. Mathematical ContentThis work makes contributions to mathematical phylogenetics in optimization, combinatorics, and statistics. We show that any tree maximizing quartet support (the criterion underlying ASTRAL) is a refinement of the networks tree of blobs under the coalescent model. Second, we give a concise proof that whether a network has a cut-edge corresponding to a given split is determined by information in certain subcollections of its 4-taxon subnetworks (quarnets). Finally, we propose valid statistical approaches for combining p-values across multiple quarnet hypothesis tests, proving that their use with specific decreasing test levels leads to statistically consistent inference as the number of loci grows. MSC codes05C90, 60J95, 62-04, 62F07, 92D15

With the increasing availability of large scale empirical fitness landscape data, there is a need for simple yet informative null models that can be used to interpret metrics of landscape ruggedness and navigability. A natural choice of a null model that maximizes ruggedness in a statistical sense assigns independent and identically distributed fitness values to the genotypes, a setting often referred to as the House-of-Cards (HoC) or mutational landscape model. In this work we examine the navigability of these landscapes, as quantified by the mean size of the adaptive basins of local fitness peaks. The adaptive basin is the set of genotypes from which a peak can be reached via selectively accessible, i.e., strictly fitnessincreasing mutational paths. Building on recent rigorous results on the statistics of accessible paths, we show that the adaptive basins in the HoC landscape encompass a positive fraction of all genotypes that is an analytically computable, increasing function of the number of alleles per site. For the four letter nucleotide alphabet, an average peak basin contains 52.8 % of all genotypes. When conditioned on peak fitness, the expected basin size increases linearly with fitness rank. The exact results on adaptive basins are complemented by an approximate analysis of gradient basins formed by greedy adaptive paths which maximize the fitness increase in each step. We argue that recent reports of large adaptive basins in empirical fitness landscapes should be reinterpreted in the light of our findings.

The relationship between genotypic and phenotypic variation is determined by the complex interaction of genetic and environmental factors. While statistical methods capable of detecting such interactions exist, an axiomatic mathematical framework that seamlessly describes the combined effects of genetic modifications and environmental exposures on a common scale is lacking. In this report, buffering concepts are used to construct a measurement system that enables the geometric representation of both gene-by-gene and gene-by-environment interactions on the extended complex plane (i.e., as projections on the Riemann sphere). In this manner, any such interaction, or combination thereof, can be precisely defined and quantified as the deviation from the neutral value calculated through the applicable complex transformation. When thus conceptualized, the framework's parameterization defines the "state space" of a given measurable phenotype along both the real and imaginary dimensions, thus establishing an unambiguous and broadly applicable method for determining the phenotypic value expected upon combinatorial changes in genetic and/or environmental variables. Remarkably, by applying these methods, it is possible to quantify the effects of any gene-by-environment interaction using the equation, AGxE=Im([z]obs*zexp)/2, where zobs and zexp are complex numbers representing the observed and expected phenotypes of a given genotype expressed in terms of the buffering parameters, and b.

Minimum viable populations (MVPs) are population levels large enough to surmount risk from demographic, environmental, and genetic stochasticity. MVPs are estimated by biologists to guide conservation practices. However, MVPs are generally estimated for a target population without regard for how they interact with intra- and inter-species population dynamics in the broader ecological community. Thus, how and why population dynamics interact with MVPs imposed by conservation biologists remain unclear. When MVPs are imposed on a continuous population model, traditional analyses fail to capture the range of possible outcomes those MVPs create. Here, we describe viability space decomposition (VSD) as a mathematical tool to systematically analyze the potential crossing of MVPs during population dynamics. We demonstrate that different extinction and survival outcomes can be recovered from a model with imposed MVPs using three VSD concepts in junction with a traditional phase portrait: mortality manifolds which separate conditions that lead to different existential outcomes, ordering manifolds which determine the order of extinction events for multiple populations, and collapse manifolds which determine the survival or extinction of one species given the loss of another. We employ these methods with a standard consumer-resource model, and the methods can be scaled to systems with more species. VSD is a useful tool for conservation biologists and community ecologists concerned with boundary crossing problems in any dynamical system.

Eusocial Hymenoptera present the highest known recombination rates among metazoans, which evolved several times independently among bees, ants and wasps. Several hypotheses have been proposed to explain this observation, including stronger selection for recombination caused by coevolving parasites and pathogens, and strong sexual selection among haploid males due to male-biased sex ratios among reproductive individuals. In this article, we explore the effects of haplodiploidy and differential selection between sexes on the evolution of recombination, by analyzing a three-locus model in which selection for recombination stems from negative epistasis between selected loci. Our analytical predictions are compared with the results of individual-based simulations in which deleterious mutations occur along a linear chromosome. Our results show that, at mutation-selection balance for deleterious alleles, increasing the strength of selection against deleterious alleles (due to the effect of male haploidy and/or sexual selection) tends to reduce selection for recombination. However, an increase in the overall magnitude of negative epistasis (which may also be due to male haploidy and/or sexual selection) combined with the fact that recombination only occurs in females may increase selection for recombination substantially. Our model also shows that, in conditions favoring recombination, increasing recombination in meioses leading to parthenogenetic ovules (and male offspring) may yield stronger benefits than in meioses leading to fertilized ovules (and female offspring).

An outbreak of New World screwworm has recently been spreading across Mexico, after more than 30 years of absence. The sterile insect technique, which consists of the massive release of sterilized males, has proven to be one of the most efficient methods for controlling the screwworm pest. However, given the limited number of sterile males available, improving the release strategy is critical. We propose a mathematical model of population dynamics adapted to the biology of Cochliomyia hominivorax and derive a feedback control function to determine the number of sterile males to release. We further construct a Luenberger observer to estimate wild fly populations from infected animal counts--the variable monitored by Mexican sanitary authorities--enabling field implementation of the control function. We show that eradication is achievable within approximately 60-100 weeks and that eradication time is governed primarily by the intrinsic biology of the system rather than by infestation magnitude. We then extend the model to a spatially explicit framework and show that when sterile male releases are applied at the outbreak focus and within a 120 km radius, eradication of the pest is attainable.

The basic reproduction number R0 is central to malaria epidemiology, yet it is typically treated as a static quantity derived under memoryless assumptions for mosquito demography. In natural systems, however, mosquito populations are shaped by delayed processes such as larval development and density-dependent feedback, introducing biological memory into vector dynamics. We develop a minimal delay-based framework that incorporates this memory into the Ross-Macdonald model by describing adult mosquito abundance with a retarded differential equation. This formulation induces a time-dependent transmission potential R0(t). Using complex analysis and the argument principle, we derive an explicit stability threshold [Formula], which separates stable from oscillatory transmission regimes. Near this threshold, delayed feedback produces slow relaxation times and sustained transient oscillations, implying that transmission potential may vary intrinsically even in the absence of external forcing. To account for ecological variability, we extend this deterministic condition into a probabilistic framework and define the stability probability as [Formula]. Numerical simulations and global sensitivity analysis show that recruitment and developmental delays are the primary drivers of instability, while adult mortality has a weaker stabilizing effect. These results indicate that malaria interventions may influence not only the magnitude of malaria transmission but also its dynamical stability. By linking delay dynamics, transmission theory, and uncertainty quantification, this framework provides a basis for stability-aware modeling and interpretation of malaria transmission under ecological variability. Author summaryMalaria transmission is often summarized by a single number, R0, treated as a fixed indicator of whether transmission will increase or decline. This assumes mosquito populations respond instantly to environmental conditions. In reality, mosquitoes develop through stages where larval conditions, such as crowding, nutrition, or temperature, affect adult populations only after a delay. This creates biological memory: todays mosquitoes reflect past environments. We show that this memory can fundamentally reshape transmission dynamics. When developmental delays are included, transmission potential is no longer constant but can fluctuate over time, even in stable environments. These fluctuations can persist or amplify depending on the balance between mosquito growth, mortality, and delay. As a result, variability in mosquito abundance or malaria transmission may arise from intrinsic dynamics rather than external drivers alone. Under ecological variability, stability becomes probabilistic, allowing estimation of how likely transmission is to remain stable. Interventions that reduce larval productivity or increase adult mortality may therefore both lower transmission and make it more predictable, improving interpretation and control strategies.

The ecological and evolutionary dynamics of populations, including viral populations, are known to be jointly shaped by deterministic and stochastic processes. While the impact of stochastic processes has been rigorously explored for viral dynamics at the level of the host population, most dynamic models for acutely-infecting respiratory viral pathogens at the within-host scale remain deterministic in their formulation. While this may be reasonable for identifying key processes shaping their within-host viral population dynamics, recent studies indicate that stochastic processes need to be invoked for understanding patterns of within-host viral evolution. Specifically, several studies have shown that viral allele frequencies can change dramatically over the time course of days in acute infections. Here, we use stochastic dynamic models to explore the role of environmental noise in shaping observed patterns of virus evolution in acute respiratory virus infections. We summarize ways in which environmental stochasticity can be biologically realized in these acute viral infections and describe within-host models that can be implemented to jointly yield viral population dynamics and evolutionary dynamics. We further develop a statistical approach to estimate the extent of environmental noise from observed within-host allele frequency changes. We test this approach on simulated data and apply it to existing influenza A virus and SARS-CoV-2 within-host data. With these applications, we show that environmental stochasticity can parsimoniously reproduce key features of empirically observed allele frequency changes without needing to invoke demographic stochasticity or to adopt Wright-Fisher model formulations with a constant effective population size. Finally, we show that purifying selection and positive selection can both still contribute to within-host viral evolution in the context of a noisy environment, providing theoretical support for studies that have found purifying and positive selection in acutely-infecting respiratory virus populations.

Gene drives, genetic constructs that can spread deleterious alleles in wild populations, have the potential to address some of the major pressing challenges of the Anthropocene such as invasive species, spread of disease vectors, and agricultural pests. However, responsible and effective deployment of gene drive requires taking into account the complex nature of real-world population connectivity networks. In particular, it is unclear how the topological position of the deployment site affects the spread process and its final outcome. Here we develop a framework for modeling gene drive spread in population connectivity networks, and study the eco-evolutionary dynamics of gene drive spread under complex population structures. We investigated the relationship between the position of the deployment site in the topology of the network and whether the gene drive is eventually lost, fixed, or maintained at an intermediate frequency. We identified network centrality measures of deployment sites that are highly correlated with the outcome of deployment for different gene drive designs and across diverse network topologies. We also show that there is a trade-off between the time-to-fixation and the final outcome, implying that multiple centrality measures of the deployment site would need to be considered when aiming to achieve rapid and successful population control using gene drives.

Infectious disease time series often show signs of epidemic transitions, such as the peaks and troughs of the time series. In these time series, key system parameters can lead to catastrophic changes in the dynamical system behaviour (often called critical transitions). Modellers have increasingly shown that early warning signals can anticipate these transitions, both critical and non-critical, in infectious disease time series. Existing methods, however, generally focus on univariate time series data, or ignore spatiotemporal patterns that may be present as a disease spreads through a population. Recent ecological literature developments expand existing temporal and spatial methods to consider the covariance matrix of multiple, related time series. However, many of these proposed signals still make an assumption of stationary time series/system equilibrium. Whilst often true in ecological modelling, disease systems are seldom at equilibrium. In this paper, we propose the usage of the eigendecomposition of the non-stationary covariance matrix as a more suitable early warning signal for epidemiological data. We first analyse the expected trends in the eigenvalues and eigenbasis of the covariance matrix on approach to a transition. Next we apply these methods to a spatially-structured susceptible-infectious-recovered model to explore how the eigenbasis may provide extra information to modellers. Finally, we test these methods on SARS-CoV-2 case data during the 2020-2021 pandemic period in England.

The general random walk model (GRW) of Hunt (2006) is used to infer directional evolution in mean trait values from sparse fossil data by modeling phenotypic change as the accumulated result of small steps with mean step sizes and step variances. Using simulations and real data cases, Ergon (2026) showed that the step variances can be estimated reasonably well only when the mean trait values have small measurement errors, while for fossil data with realistic measurement errors they appear to be extremely difficult to find, and they are often found to be negative. In the simulations Ergon (2026) assumed that the true phenotypic mean values were known. Here, I essentially repeat these simulations under the assumption that only mean trait values with large measurement errors are known, and based on weighted mean squared error (WMSE) comparisons the conclusion is that weighted least squares (WLS) is a better method than GRW. A second conclusion is that WLS is a better method also in the possibly rare cases with large measurement errors where the GRW parameters are estimated well. The GRW method is simply not flexible enough to handle such cases. A third conclusion is that Akaike Information Criterion (AIC) results for GRW models with large measurement errors relative to the step variance may be overly optimistic.

Multicellular organisms regularly encounter microbes, which are, however, only rarely pathogenic. Our understanding of this phenomenon is currently restricted due to lacking theory on evolutionary transitions between non-pathogenic and pathogenic microbial lifestyles. Here we addressed this gap by investigating a mathematical model of host-microbe interactions that is based on the danger theory of immunology, which states that danger signals related to host tissue damage play a key role in activating immune responses. We formally implemented this idea by assuming that immune activation increases with costs that microbes cause to their host, and we compared this to scenarios in which immune activation depends only on the presence or load of infecting microbes. Our model analysis revealed that cost-based - but not presence or load-based - immune activation favours the evolution of avirulence and associated non-pathogenic microbial lifestyles. Based on our results, we propose the danger hypothesis of virulence evolution which states that evolution towards avirulence and intermediate virulence are both possible - depending on whether hosts can accurately assess costs generated by microbes. The idea that basic host immune responses can select for avirulence offers a new explanation for why most microbes are not pathogenic to a given host.

We describe a statistical agent-based model (SABM) for binary spatio-temporal data in which the occupancy of each cell evolves as a Bernoulli mixture of three mechanistically distinct processes: local persistence, anisotropic neighborhood dispersal, and long-distance dispersal. The model is embedded in a hierarchical Bayesian framework with conjugate Beta full-conditionals for the persistence and long-distance parameters and a Dirichlet prior on the directional dispersal kernel. A nonstationary extension links the dispersal kernel to a latent habitat-suitability surface through directional gradients of a Gaussian process. We show that, in the small-step regime, the Lagrangian recurrence for the dispersal kernel scales to a classical two-dimensional advection-diffusion partial differential equation whose drift and dispersion coefficients are the first and second moments of the dispersal probabilities. We provide an MCMC algorithm exploiting the exact full-conditionals and demonstrate parameter recovery and PDE-scaling agreement in a simulated example.